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ABSTRACT 

Aims. To investigate the stability and angular momentum transport by the strato-rotational instability in the nonlinear regime. 
Methods. The hydrodynamic compressible equations are solved in a cartesian box in which the outer cylinder is embedded. 
Gravity along the rotation axis leads to density stratification. No-slip boundary conditions are used in the radial direction, 
while free-slip conditions are used on the two ends of the cylinders. 

Results. The strato-rotational instability is confirmed and the Reynolds stress is shown to transport angular momentum away 
from the axis. However, the growth rate decreases with increasing Reynolds number. This, as well as the presence of boundaries 
renders this instability less relevant for astrophysical applications. 
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1. Introduction 

The magneto-rotational instability (MRI) is now com¬ 
monly invoked in order to explain the origin of turbulence 
in accretion discs (Balbus & Hawley 1998). This is because 
there exists no purely hydrodynamic local instabilities in 
accretion discs in the sense that periodic shearing box cal¬ 
culations do not yield instability. Already when vertical 
shear is included, there is a linear instability, but for thin 
discs its growth rate is small compared with that of the 
MRI (Urpin & Brandenburg 1998; Arlt & Urpin 2004). 
Purely hydrodynamic instabilities are of interest for proto- 
stellar discs where the conductivity is usually so poor that 
magnetic effects are unimportant, and hence the MRI may 
be irrelevant. Another possibility for a hydrodynamic in¬ 
stability is the nonlinear shear instability (Chagelishvili et 
al. 2003, Afshordi, Mukhopadhyay & Narayan 2005). The 
general difficulty with such models has been discussed by 
Balbus, Hawley & Stone (1996). However, a conclusive res¬ 
olution of this problem using simulations remains difficult 
because of the large Reynolds numbers required. 

In an attempt to address the possibility of hydrody¬ 
namic instabilities in accretion discs, one can resort to ex¬ 
periments. All these experiments have boundaries, so their 
relevance to accretion discs is questionable. Nevertheless, 
such experiments remain interesting in their own right. 
The possibility of a nonlinear instability of a Taylor- 
Couette flow experiments has recently been addressed by 


Richard & Zahn (1999). The issue was revived by re¬ 
cent findings of Molemaker, McWilliams & Yavneh (2001) 
who found the possibility of a linear instability in Taylor- 
Couette flow with density stratification along the axis. 
Dubrulle et al. (2005) confirmed this instability for a range 
of different radial boundaries conditions. 

More recent work by Shalybkov & Rudiger (2005), 
who considered regular Taylor-Couette flow in the pres¬ 
ence of no-slip boundaries, gave a detailed stability cri¬ 
terion for the instability, which is henceforth referred 
to as the strato-rotational instability (SRI). Withjack 
& Chen (1974) observed nonaxisymmetric patterns in 
density-stratified Taylor-Couette flows close and beyond 
the Rayleigh line. 

Umurhan (2005) shows that, within the framework of 
the quasi-hydrostatic semi-geostrophic approximation, the 
instability only survives in the presence of no-slip radial 
boundaries. The purpose of the present paper is to demon¬ 
strate the existence of this instability using direct simula¬ 
tions and to calculate its saturation level in the nonlinear 
regime. In contrast to nonlinear instabilities, the SRI can 
be studied conveniently by simulations, because the crit¬ 
ical Reynolds numbers are well below 10^. We are thus 
also able to compute sign and, in the nonlinear regime, 
also the amount of angular momentum transport due to 
the SRI. Let us emphasize here that the sign of angular 
momentum transport (inwards or outwards) is nontrivial. 
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There are several examples in the literature (mostly for 
rotating convection) where the cross correlation 

Qrtj) — (^) 

of the one-point correlation tensor 

Qzj = ( 2 ) 

turns out to be negative, i.e. angular momentum is trans¬ 
ported inwards by the rotating turbulence (see Ryu & 
Goodman 1992; Kley, Papaloizou & Lin 1993; Stone & 
Balbus 1996). The primes in Q and © denote the de¬ 
partures from the toroidally averaged flow, i.e. u' = u — u, 
where 

u = u d(j> / 2 tt . (3) 

Here we shall find that the zonal flux of angular momen¬ 
tum, o, always proves to be positive i.e. the angular mo¬ 
mentum is transported outwards by the SRI. 


where ^in(T’) = 1 when |r| < i?in and 0 otherwise, 
and Cout(»’) = 1 when |r| > i?out and 0 otherwise, r 
is a relaxation time (comparable to the time step), and 
Min = (J^in — J^out) X r is the linear velocity of the inner 
cylinder in the rotating frame. The temperature is given 
by CpT = Cs /(7 — 1). Gravity is written as g = (0,0, —g), 
where the value of g is varied. 

The initial stratihcation is isothermal (cg = Cso = 
const) with Inp = —z/iL and s/Cp = (7 — \)zlH^ where 
H = Cso/( 75 ) is the density scale height. The Brunt- 
Vaisala frequency N is given by 


= 


Cp dz 


= (7-l)^= 1 - 


9_ 

H’ 


( 11 ) 


which, following Shalybkov & Rudiger (2005), is expressed 
in terms of the Froude number, Fr = flin/N. For 7 = 5/3 
and g = 1 we have N = 0.27 flin, so Fr = 3.7. 

Viscosity and thermal diffusivity are specified in terms 
of Reynolds and Prandtl numbers, 


2. The model 


Re = Hi„i?in(Rout - Rin)/j^, Pr = i^/x, (12) 


We consider the governing equations for a compressible 
perfect gas. 


Dlnp 

Dt 


— V • u, 


( 4 ) 


Dm 

~Dt 


= g - c^{V \np + Vs/Cp) + 


. + F, 


drive 


( 5 ) 


respectively. In all cases considered below we assume Pr = 
10 , i.e. thermal diffusion is small compared with viscosity. 
(We need to keep a certain minimum amount of thermal 
diffusivity to stabilize the code.) 

We assume the radial boundaries to be insulating, i.e. 
dT/dm = 0 , and for the velocity we adopt a no-slip bound¬ 
ary condition, i.e. 


= V ■ (cpPxVT) , ( 6 ) 

where u is the velocity, p the density, s the specific entropy, 
g = const is gravity, Cp and c„ are the specific heats at con¬ 
stant pressure and volume, respectively, 7 = Cp/c„ = 5/3 
is their assumed ratio, Cg is the sound speed, T is the tem¬ 
perature, and the two are related to the other quantities 
via 

(7 - l)cpT = cl = c^o exp [ 7 s/cp -b (7 - 1) In p/po]. (7) 

Here, Cgo and po are normalization constants. These con¬ 
stants also define the zero point of the specific entropy. 
The viscous force is 

T’visc = (V^M-b 4VV • M-f 2S • V Inp) , ( 8 ) 

where = 4(uij-bMj,^) — V-m is the traceless rate of 
strain tensor, where commas denote partial differentiation. 
We consider a rotating frame of reference, such that in 
that frame the angular velocity of the outer cylinder is 
zero. The resulting apparent force is 

FYot — 2i^out X M -b (9) 

where m = (x, y, 0) is the cylindrical radius vector. In our 
frame of reference, the angular velocities of the outer and 
inner cylinders are imposed by a body force of the form 

T'drive — 'C [(m '?^in)^in(^) “b M^out(^)] , (4^) 


M = Min (on |r| = i?in), M = 0 (on |r| = Rout)- (13) 

In the vertical direction we adopt insulating free-slip 
boundary conditions. The vertical extent of the cylindrical 
shell is in all cases = 4. It is convenient to normalize 
the various quantities with respect to their values on the 
inner cylinder and to introduce the ratios 

g = Rin/Rout, A = Hout/Hin- (14) 

Following Shalybkov & Rudiger (2005), the flow is strato- 
rotationally unstable when 

ff' < pL < fj (SRI, nonaxisymmetric). (15) 

For p < if the flow is always Rayleigh unstable both 
to axisymmetric and to nonaxisymmetric perturbations. 
In the following we consider the cases fj = 0.5 and 0.25 
(Rout = 2 and 4, using Ri^ = 1), so we expect the SRI to 
occur in the ranges 0.25 < p < 0.5 and 0.06 < p < 0.25 in 
the two cases. For p < 0.25 (or 0.06) the flow is Rayleigh- 
unstable. 

For the runs presented below the maximum flow speed 
is at the most 1.2 in units where Din = Rin = 1. In order 
that the flow remains everywhere subsonic, i.e. max(|M|) < 
Cg, we choose Cgo = 3, which means that the sound speed 
varies around 3. For the SRI compressibility is however 
not essential. 

We model the problem by embedding the cylinder 
into a box where the velocities corresponding to rigid 
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Fig. 1. Amplitudes for /t = 0, 0.2, 0.3, and 0.6. Here, 
g/(i?inf2fn) = 1 = 500. 



Fig. 2. Amplitudes for gj = 1, 2, 5, and 10, using 
p, = 0.3 and Re = 500. 


inner and outer cylinders is prescribed (see Dobler, 
Shukurov & Brandenburg 2002 for details regarding this 
approach). We solve the equations in cartesian geometry 
using the Pencil Code^, which is a memory-efficient sixth- 
order finite difference code using the 2WRK3 scheme of 
Williamson (1980). For most of the cases we consider a 
Reynolds number of 500 where a resolution of 128^ mesh- 
points is sufficient. In a few cases we considered larger 
Reynolds numbers where a higher resolution of up to 512^ 
meshpoints was used. 

3. Results 

3.1. Dependence on fi 

As a simple diagnostics of the instability we use the root 
mean square velocity, Uz^rms- In Fig. ^ we plot the evo¬ 
lution of Uz,rms for different vafues of p. As expected, 
no instabifity is found for p = 0.6. For smaffer vafues 
(p = 0 and 0.2), the ffow is Rayfeigh unstable. For p = 0.3 
the flow is strato-rotationally unstable (Fig. Since 
such flows are stable to axisymmetric perturbations by 
the Solberg-Hpiland criteria (Rudiger, Arlt & Shalybkov 
2002), the instability must be necessarily nonaxisymmet- 
ric, which is indeed the case^. 

3.2. Dependence on g 

For weak stratification (small values of g) the growth rate 
of the instability is small. For small gap width, Re = 500 
and g/RinD.f^ = 2 we find the largest rms velocity. For 
= 5 the growth rate is maximum while for g = 10 
no instability is found; see Fig.|21 However, when Re is in¬ 
creased to 1000, also the case with g = 10 yields instabil¬ 
ity. The number of azimuthal cells increases as the value 
of g is increased; see Fig. 01 where we show the velocity 

^ http://www.nordita.dk/software/pencil-code 

^ Animations of the flow, showing also the nonax- 
isymmetric nature of the instability, can be found on 
http://www.nordita.dk/'brandenb/movies/couette/ 


field in a meridional cross-section. Note that for g = 5 
there are 4 eddies in the z direction. Since Lz/Rm = 4, 
the vertical wavelength is 1. 

3.3. Critical Reynolds number 

The growth rates are on the order of (0.01...0.02) x llin, 
which is relatively small. The dependence of the growth 
rate A on the Reynolds number Re is shown in Fig. 0) 
For Re < 400 the vertical rms velocity seems to grow in 
an oscillatory fashion. However, for Re > 500, the growth 
rate decreases with increasing Reynolds number. 

3.4. Reynolds stress 

In order to estimate the ability of the flow to transport 
angular momentum, we compute the radial component 
of the Reynolds stress. It turns out that Qnj) is positive 
throughout most of the domain, so the flow is transport¬ 
ing angular momentum outwards. The normalized stress, 
Qr^/iRin^in) is around 5 x 10“"^ for g/{RinD,f^) around 
2 and 5. For smaller values of g, the vertical scale of the 
eddies becomes comparable to or exceeds the vertical ex¬ 
tent of the cylinders, so the SRI begins to be suppressed. 
Likewise, for larger values of g, the vertical scale of the 
eddies becomes comparable to the scale where viscosity 
effects can still be regarded as weak. 

The fact that the ratio Qr,p/{RfnDf^) is much less than 
unity suggests that angular momentum transport is rela¬ 
tively inefficient. This supplements the earlier impression 
based on the small value of the normalized growth rate of 
the instability. 

Let us now express the Reynolds stress in astrophys- 
ically relevant terms by normalizing it with respect to 
where we have chosen the inner radius as ref¬ 
erence. (We note that in accretion discs, iLinllin is pro¬ 
portional to the sound speed.) The Reynolds stress nor- 



















4 


A. Brandenburg and G. Rudiger: The angular momentum transport by the strato-rotational instability 





Re 

Fig. 4. Dependence of the growth rate A on the Reynolds 
number Re. The critical value is about Re « 400. For 
Re > 500, A decreases with increasing Re. 

Table 1. Dependence of the Reynolds stress on gravity 
for container with small gap. 


g/{Rin^in) 

Q.^/(10-"Rf„DfJ 

g../(io-"Rf„Drj 

0.5 

> 0.16 

> 0.04 

1 

0.9 

4.6 

2 

4.3 

5.4 

5 

4.2 

1.1 


Table 2. Dependence of the Reynolds stress on grav¬ 
ity for the container with large gap. Note that for 
g/{Rini^fn) — 2 the perturbations are slowly decaying. 


9 / {R-A^) = 5 



- 2-10 1 2 
X / Rin 


Fig. 3. Velocity pattern in the meridional plane for dif¬ 
ferent values of and Re = 500. 


ff/(RinDrj g../(io-^ii;rju 

2 0.1 0.02 

5 3 0.6 


Using the data from Tabled we find 

ass = (1-5) X 10-4 . (17) 

One may be tempted to insert here the standard result for 
Shakura-Sunyaev discs, H/R = 0.03, which would suggest 
rather large values of the stress normalized in this way; 
ass = 0.1...0.5. However, this argument assumes that the 
instability survives for small values of H/R, which is un¬ 
fortunately not true. Furthermore, for larger radial the 
instability seems to develop only near the inner cylinder; 
see Fig. 121 Nevertheless, provided the flow is unstable, the 
resulting normalized stress, Qr^/{Rfn^fn)^ is still roughly 
the same. 


malized in this way is referred to as the Shakura-Sunyaev ^ Conclusions 
alpha. 


ass = Qr<^/ (77in^?n) • 


In the present work we have explored the nonlinear sat- 
(16) uration of the strato-rotational instability using a finite 
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Fig. 5. Velocity pattern in the meridional plane for dif¬ 
ferent values of (For < 2 the per¬ 

turbations are slowly decaying.) 


difference method. We verify the stability criterion of the 
Taylor-Couette flow in the presence of vertical stratifica¬ 
tion. In agreement with earlier results by Shalybkov & 
Rudiger (2005), a nonaxisymmetric instability is found 
when the azimuthal velocity on the inner cylinder, 
exceeds that on the outer cylinder, i?outf^out- This leads to 
instability in a regime that would be stable by the Solberg- 
Hpiland criteria, which only apply to axisymmetric flows 
on in the tight-winding approximation (Rudiger, Arlt & 
Shalybkov 2002). 

Increasing Re sufficiently beyond the critical value for 
the instability, the growth rate decreases with increasing 
Re. This is related to the fact that the flow is driven via 
the viscous force, which is obviously not relevant in as¬ 
trophysics, but it is relevant to the experiments that are 
currently being designed to address these questions that 
are astrophysically motivated in other ways. 
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